### 抓所需的資料
ticker = '^GSPC'
start = '2017-01-01'
end = as.Date(Sys.time())
# S&P500 報酬率
# Asset = getSymbols(ticker , src = 'yahoo' , auto.assign = FALSE , from = start , to = end)
# save(Asset , file = 'C:/Users/amyhs/Desktop/碩士課程/時間序列/期中/Asset.RData')
load('C:/Users/amyhs/Desktop/碩士課程/時間序列/期中/Asset.RData')
Asset$Asset.LogReturn = dailyReturn(Ad(Asset) , type = 'log')
# 加入星期一為解釋變數
Asset$Asset.Mon = NA
Asset$Asset.Mon = xts(x = as.integer(weekdays(index(Asset)) == '星期一') , order.by = index(Asset))
# 若有兩個解釋變數用cbind合併 再轉為matrix格式(regressor規定用matrx)
# mean:直接帶入;var:帶入要加平方
Regressor_mean = as.matrix(Asset$Asset.Mon)
Regressor_var = as.matrix(Asset$Asset.Mon)
因後面題目要考慮風險溢酬及波動不對稱性,故此題先加入in-mean及gjrGARCH效應,但實際上並不會如此,因為容易造成模型overfitting。 建立ARMA-gjrGARCH-in-mean模型的步驟為先決定ARMA的order,再決定gjrGARCH-in-mean的order。
Asset_2017 = Asset['2017']
Return_2017 = Asset_2017$Asset.LogReturn
plot(Return_2017 , type = 'l' , main = ticker , ylab = 'Log Return')
由圖判斷資料為定態。
par(mfrow = c(2,1) , mai=c(0.5,0.5,0.7,0.5))
acf(Return_2017 , lag.max = 20 , main = sprintf("%s - Log Return 2017 (ACF)" , ticker))
pacf(Return_2017 , lag.max = 20 , main = sprintf("%s - Log Return 2017 (PACF)" , ticker))
由圖可發現ACF、PACF皆為tail-off,故假設模型為ARMA(0,0)
par(mfrow = c(2,1) , mai=c(0.5,0.5,0.7,0.5))
acf(Return_2017^2 , lag.max = 20 , main = sprintf("%s - Square of Log Return (ACF)" , ticker))
pacf(Return_2017^2 , lag.max = 20 , main = sprintf("%s - Square of Log Return (PACF)" , ticker))
將報酬率平方後可以看到ACF為cut-off,PACF為tail-off,表示此時間序列有Garch效應
DataSize = length(Asset$Asset.LogReturn) # 資料大小
WindowSize = length(Return_2017)
OutSample = DataSize - WindowSize
ARMA_p = 0 # AR order
ARMA_q = 0 # MA order
spec1 = arfimaspec(mean.model = list(armaOrder = c(ARMA_p, ARMA_q),
include.mean = TRUE),
distribution.model = "norm")
# 加速收斂
solver_control = list(tol=1e-5, delta=1e-5, trace=1)
modelfit1 = arfimafit(spec = spec1,
data = Return_2017,
solver = "hybrid",
solver.control = solver_control
)
##
## Iter: 1 fn: -1018.9354 Pars: 0.000670 0.004176
## Iter: 2 fn: -1018.9354 Pars: 0.0006701 0.0041757
## solnp--> Completed in 2 iterations
modelfit1
##
## *----------------------------------*
## * ARFIMA Model Fit *
## *----------------------------------*
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000670 0.000264 2.5423 0.011012
## sigma 0.004176 0.000186 22.4060 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000670 0.000214 3.1304 0.001746
## sigma 0.004176 0.000306 13.6576 0.000000
##
## LogLikelihood : 1018.935
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.1031
## Bayes -8.0750
## Shibata -8.1032
## Hannan-Quinn -8.0918
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 4.430 0.03532
## Lag[2*(p+q)+(p+q)-1][2] 4.824 0.04532
## Lag[4*(p+q)+(p+q)-1][5] 5.471 0.11947
##
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1443 0.7041
## Lag[2*(p+q)+(p+q)-1][2] 0.1796 0.8669
## Lag[4*(p+q)+(p+q)-1][5] 1.0589 0.8462
##
##
## ARCH LM Tests
## ------------------------------------
## Statistic DoF P-Value
## ARCH Lag[2] 0.2024 2 0.9037
## ARCH Lag[5] 3.7154 5 0.5911
## ARCH Lag[10] 6.0199 10 0.8136
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 0.1261
## Individual Statistics:
## mu 0.03802
## sigma 0.10359
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 0.61 0.749 1.07
## Individual Statistic: 0.35 0.47 0.75
##
##
## Elapsed time : 0.01736712
modelfit1_residual = modelfit1@fit$residuals
par(mfrow = c(2,1) , mai=c(0.5,0.5,0.7,0.5))
acf(modelfit1_residual , main = "modelfit1 - residual")
pacf(modelfit1_residual , main = "modelfit1 - residual")
Weighted.Box.test(modelfit1_residual , lag = 10, type="Ljung-Box")
Weighted.Box.test(modelfit1_residual , lag = 20, type="Ljung-Box")
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit1_residual
## Weighted X-squared on Residuals for fitted ARMA process = 7.0888,
## Shape = 3.9286, Scale = 1.4000, p-value = 0.2449
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit1_residual
## Weighted X-squared on Residuals for fitted ARMA process = 11.885,
## Shape = 7.6829, Scale = 1.3667, p-value = 0.3192
圖中可看出模型ARMA(0,0)殘差的ACF、PACF皆顯示無相關,且Weighted Ljung-Box test通過,故判斷殘差為white noise。
| 由以上過程決定使用ARMA(0,0)模型。 |
決定gjr-GARCH的order(窮舉法)
ARMA_p = 0
ARMA_q = 0
GARCH_p = 0 # GARCH order
GARCH_q = 0 # ARCH order
spec2_1 = ugarchspec(mean.model = list(armaOrder = c(ARMA_p , ARMA_q),
include.mean = TRUE,
archm = TRUE , # 是否要做in-mean Model
external.regressors = NULL),
variance.model = list(model = "gjrGARCH" ,
garchOrder = c(GARCH_q , GARCH_p),
variance.targeting = TRUE , # 較易收斂
external.regressors = NULL) ,
distribution.model = "norm")
solver_control = list(tol=1e-5, delta=1e-5, trace=1)
modelfit2_1 = ugarchfit(spec = spec2_1,
data = Return_2017,
solver = "hybrid",
solver.control = solver_control
)
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
##
## Iter: 1 fn: 1.1000 Pars: 0.0006734 NA
## solnp--> Solution not reliable....Problem Inverting Hessian.
##
## Trying nlminb solver...
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
## 0: 1.1000000: 0.000673408 nan
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
##
## Trying gosolnp solver...
##
## Calculating Random Initialization Parameters...ok!
##
## Excluding Inequality Violations...
##
## ...Excluded 500/500 Random Sequences
##
## Evaluating Objective Function with Random Sampled Parameters...ok!
##
## Sorting and Choosing Best Candidates for starting Solver...ok!
##
## Starting Parameters and Starting Objective Function:
## [,1]
## par1 NA
## par2 NA
## objf NA
##
## gosolnp-->Starting Solver
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
##
## Iter: 1 fn: -914.1421 Pars: NA NA
## solnp--> Solution not reliable....Problem Inverting Hessian.
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
##
##
## arfimafit-->warning: Error in try(.C("gjrgarchfilterC", model = as.integer(modelinc[1:21]), :
## 外部函數呼叫時不能有 NA/NaN/Inf (引數 2)
modelfit2_1
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : gjrGARCH(0,0)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Convergence Problem:
## Solver Message: Error in is.nloptr(ret) : x0 contains NA
##
在fit GARCH(0,0)模型時出現收斂問題(Convergence Problem) Solver Message: Error in is.nloptr(ret) : x0 contains NA,故此模型不適合。
ARMA_p = 0
ARMA_q = 0
GARCH_p = 0 # GARCH order
GARCH_q = 1 # ARCH order
spec2_2 = ugarchspec(mean.model = list(armaOrder = c(ARMA_p , ARMA_q),
include.mean = TRUE,
archm = TRUE , # 是否要做in-mean Model
external.regressors = NULL),
variance.model = list(model = "gjrGARCH" ,
garchOrder = c(GARCH_q , GARCH_p),
variance.targeting = TRUE , # 較易收斂
external.regressors = NULL) ,
distribution.model = "norm")
solver_control = list(tol=1e-5, delta=1e-5, trace=1)
modelfit2_2 = ugarchfit(spec = spec2_2,
data = Return_2017,
solver = "hybrid",
solver.control = solver_control
)
##
## Iter: 1 fn: -1019.1396 Pars: 0.001385 -0.157271 0.032357 -0.064251
## Iter: 2 fn: -1019.1406 Pars: 0.001373 -0.155505 0.031344 -0.062292
## solnp--> Completed in 2 iterations
modelfit2_2
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : gjrGARCH(1,0)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001373 0.001100 1.24749 0.21222
## archm -0.155505 0.254736 -0.61045 0.54156
## alpha1 0.031344 0.086263 0.36335 0.71634
## gamma1 -0.062292 0.094741 -0.65750 0.51086
## omega 0.000018 NA NA NA
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001373 0.001723 0.79685 0.42554
## archm -0.155505 0.413214 -0.37633 0.70667
## alpha1 0.031344 0.066371 0.47225 0.63675
## gamma1 -0.062292 0.069691 -0.89383 0.37141
## omega 0.000018 NA NA NA
##
## LogLikelihood : 1019.141
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.0888
## Bayes -8.0326
## Shibata -8.0893
## Hannan-Quinn -8.0662
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 4.376 0.03644
## Lag[2*(p+q)+(p+q)-1][2] 4.837 0.04498
## Lag[4*(p+q)+(p+q)-1][5] 5.522 0.11628
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.03017 0.8621
## Lag[2*(p+q)+(p+q)-1][2] 0.05627 0.9515
## Lag[4*(p+q)+(p+q)-1][5] 0.96481 0.8676
## d.o.f=1
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[2] 0.05139 0.500 2.000 0.8207
## ARCH Lag[4] 0.23269 1.397 1.611 0.9502
## ARCH Lag[6] 2.08525 2.222 1.500 0.6612
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 0.4968
## Individual Statistics:
## mu 0.01626
## archm 0.04738
## alpha1 0.34808
## gamma1 0.22700
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.07 1.24 1.6
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.67859 0.4980
## Negative Sign Bias 0.07271 0.9421
## Positive Sign Bias 0.40004 0.6895
## Joint Effect 1.71114 0.6345
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 56.81 1.223e-05
## 2 30 72.55 1.349e-05
## 3 40 80.95 9.178e-05
## 4 50 89.44 3.714e-04
##
##
## Elapsed time : 0.8717721
Residual Analysis
modelfit2_2_std_residual = modelfit2_2@fit$residuals/modelfit2_2@fit$sigma
par(mfcol = c(2,2) , mai=c(0.75,0.75,0.7,0.7))
acf(modelfit2_2_std_residual , lag.max = 20 , main = sprintf('%s - Standardized Residual' , ticker))
pacf(modelfit2_2_std_residual , lag.max = 20 , main = '')
acf(modelfit2_2_std_residual^2 , lag.max = 20 , main = sprintf('%s - Square of Standardized Residual' , ticker))
pacf(modelfit2_2_std_residual^2 , lag.max = 20 , main = '')
Weighted.Box.test(modelfit2_2_std_residual , lag = 10, type="Ljung-Box")
Weighted.Box.test(modelfit2_2_std_residual , lag = 20, type="Ljung-Box")
Weighted.Box.test(modelfit2_2_std_residual^2 , lag = 10, type="Ljung-Box")
Weighted.Box.test(modelfit2_2_std_residual^2 , lag = 20, type="Ljung-Box")
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit2_2_std_residual
## Weighted X-squared on Residuals for fitted ARMA process = 7.1563,
## Shape = 3.9286, Scale = 1.4000, p-value = 0.2385
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit2_2_std_residual
## Weighted X-squared on Residuals for fitted ARMA process = 11.968,
## Shape = 7.6829, Scale = 1.3667, p-value = 0.312
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit2_2_std_residual^2
## Weighted X-squared on Residuals for fitted ARMA process = 3.3015,
## Shape = 3.9286, Scale = 1.4000, p-value = 0.7755
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit2_2_std_residual^2
## Weighted X-squared on Residuals for fitted ARMA process = 6.6661,
## Shape = 7.6829, Scale = 1.3667, p-value = 0.8522
ARMA_p = 0
ARMA_q = 0
GARCH_p = 1
GARCH_q = 0
spec2_3 = ugarchspec(mean.model = list(armaOrder = c(ARMA_p , ARMA_q),
include.mean = TRUE,
archm = TRUE , # 是否要做in-mean Model
external.regressors = NULL),
variance.model = list(model = "gjrGARCH" ,
garchOrder = c(GARCH_q , GARCH_p),
variance.targeting = TRUE , # 較易收斂
external.regressors = NULL) ,
distribution.model = "norm")
modelfit2_3 = ugarchfit(spec = spec2_3,
data = Return_2017,
solver = "hybrid",
solver.control = solver_control
)
##
## Iter: 1 fn: -1018.9350 Pars: 0.0004756 0.0464003 0.9000000
## Iter: 2 fn: -1018.9350 Pars: 0.0004756 0.0464003 0.9000000
## solnp--> Completed in 2 iterations
## Warning in .makefitmodel(garchmodel = "gjrGARCH", f = .gjrgarchLLH, T = T, :
## rugarch-->warning: failed to invert hessian
modelfit2_3
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : gjrGARCH(0,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000476 NA NA NA
## archm 0.046400 NA NA NA
## beta1 0.900000 NA NA NA
## omega 0.000002 NA NA NA
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000476 NA NA NA
## archm 0.046400 NA NA NA
## beta1 0.900000 NA NA NA
## omega 0.000002 NA NA NA
##
## failed to invert hessian
## LogLikelihood : 1018.935
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.0951
## Bayes -8.0530
## Shibata -8.0954
## Hannan-Quinn -8.0781
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 4.430 0.03532
## Lag[2*(p+q)+(p+q)-1][2] 4.824 0.04532
## Lag[4*(p+q)+(p+q)-1][5] 5.471 0.11947
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1443 0.7040
## Lag[2*(p+q)+(p+q)-1][2] 0.1797 0.8668
## Lag[4*(p+q)+(p+q)-1][5] 1.0589 0.8462
## d.o.f=1
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[2] 0.06964 0.500 2.000 0.7919
## ARCH Lag[4] 0.26959 1.397 1.611 0.9397
## ARCH Lag[6] 2.01417 2.222 1.500 0.6771
## Error in t.default(grad): 引數不是矩陣
在fit GARCH(1,0)模型時出現警告訊息warning: failed to invert hessian表示資料在這個模型下共線性的影響太嚴重,故此模型不適合。
ARMA_p = 0
ARMA_q = 0
GARCH_p = 1
GARCH_q = 1
spec2_4 = ugarchspec(mean.model = list(armaOrder = c(ARMA_p , ARMA_q),
include.mean = TRUE,
archm = TRUE , # 是否要做in-mean Model
external.regressors = NULL),
variance.model = list(model = "gjrGARCH" ,
garchOrder = c(GARCH_q , GARCH_p),
variance.targeting = TRUE , # 較易收斂
external.regressors = NULL) ,
distribution.model = "norm")
modelfit2_4 = ugarchfit(spec = spec2_4,
data = Return_2017,
solver = "hybrid",
solver.control = solver_control
)
##
## Iter: 1 fn: -1018.9252 Pars: 0.00017398 0.11975019 0.00002314 0.95215001 0.00156672
## Iter: 2 fn: -1018.9272 Pars: 0.00016785 0.12122922 0.00000386 0.95154043 0.00176747
## solnp--> Completed in 2 iterations
modelfit2_4
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : gjrGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000168 0.001289 0.130202 0.89641
## archm 0.121229 0.305513 0.396805 0.69151
## alpha1 0.000004 0.001996 0.001934 0.99846
## beta1 0.951540 0.045872 20.743408 0.00000
## gamma1 0.001767 0.023718 0.074522 0.94059
## omega 0.000001 NA NA NA
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000168 0.001880 0.089298 0.92885
## archm 0.121229 0.420895 0.288027 0.77333
## alpha1 0.000004 0.013921 0.000277 0.99978
## beta1 0.951540 0.066748 14.255806 0.00000
## gamma1 0.001767 0.031032 0.056957 0.95458
## omega 0.000001 NA NA NA
##
## LogLikelihood : 1018.927
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.0791
## Bayes -8.0089
## Shibata -8.0799
## Hannan-Quinn -8.0508
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 4.446 0.03498
## Lag[2*(p+q)+(p+q)-1][2] 4.827 0.04524
## Lag[4*(p+q)+(p+q)-1][5] 5.459 0.12027
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1467 0.7017
## Lag[2*(p+q)+(p+q)-1][5] 1.0049 0.8586
## Lag[4*(p+q)+(p+q)-1][9] 2.8717 0.7800
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.1801 0.500 2.000 0.6712
## ARCH Lag[5] 2.1992 1.440 1.667 0.4289
## ARCH Lag[7] 3.1710 2.315 1.543 0.4819
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 5.2482
## Individual Statistics:
## mu 0.06202
## archm 0.04076
## alpha1 0.65039
## beta1 0.43065
## gamma1 0.14828
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.8720 0.3841
## Negative Sign Bias 0.4118 0.6808
## Positive Sign Bias 0.1465 0.8836
## Joint Effect 1.4021 0.7050
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 47.73 2.811e-04
## 2 30 70.63 2.482e-05
## 3 40 70.75 1.391e-03
## 4 50 81.87 2.240e-03
##
##
## Elapsed time : 0.5284832
Residual Analysis
modelfit2_4_std_residual = modelfit2_4@fit$residuals/modelfit2_4@fit$sigma
par(mfcol = c(2,2) , mai=c(0.75,0.75,0.7,0.7))
acf(modelfit2_4_std_residual , lag.max = 20 , main = sprintf('%s - Standardized Residual' , ticker))
pacf(modelfit2_4_std_residual , lag.max = 20 , main = '')
acf(modelfit2_4_std_residual^2 , lag.max = 20 , main = sprintf('%s - Square of Standardized Residual' , ticker))
pacf(modelfit2_4_std_residual^2 , lag.max = 20 , main = '')
Weighted.Box.test(modelfit2_4_std_residual , lag = 10, type="Ljung-Box")
Weighted.Box.test(modelfit2_4_std_residual , lag = 20, type="Ljung-Box")
Weighted.Box.test(modelfit2_4_std_residual^2 , lag = 10, type="Ljung-Box")
Weighted.Box.test(modelfit2_4_std_residual^2 , lag = 20, type="Ljung-Box")
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit2_4_std_residual
## Weighted X-squared on Residuals for fitted ARMA process = 7.0646,
## Shape = 3.9286, Scale = 1.4000, p-value = 0.2472
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit2_4_std_residual
## Weighted X-squared on Residuals for fitted ARMA process = 11.84,
## Shape = 7.6829, Scale = 1.3667, p-value = 0.3231
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit2_4_std_residual^2
## Weighted X-squared on Residuals for fitted ARMA process = 3.2041,
## Shape = 3.9286, Scale = 1.4000, p-value = 0.7902
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit2_4_std_residual^2
## Weighted X-squared on Residuals for fitted ARMA process = 6.4531,
## Shape = 7.6829, Scale = 1.3667, p-value = 0.8693
因模型一未收斂,故無法計算BIC值
infocriteria(modelfit2_2)
infocriteria(modelfit2_3)
infocriteria(modelfit2_4)
##
## Akaike -8.088769
## Bayes -8.032587
## Shibata -8.089267
## Hannan-Quinn -8.066160
##
## Akaike -8.095100
## Bayes -8.052963
## Shibata -8.095381
## Hannan-Quinn -8.078143
##
## Akaike -8.079101
## Bayes -8.008873
## Shibata -8.079874
## Hannan-Quinn -8.050839
| 模 型 | Weighted Ljung-Box test | BIC | 是否選用 |
|---|---|---|---|
| gjr-GARCH(0,0) | 未收斂 | - | X |
| gjr-GARCH(0,1) | 通過 | -8.032587 | O |
| gjr-GARCH(1,0) | 共線性影響太嚴重 | - | X |
| gjr-GARCH(1,1) | 通過 | -8.008873 | X |
| 由上表決定使用gjr-GARCH(0,1)模型。 |
因此最終模型為ARMA(0,0)-gjrGARCH(0,1)-in-mean Model: \[ Y_t - \mu_t = a_t \] \[ a_t \mid \mathcal{F_{t-1}} ~ \stackrel{iid}{\backsim} N(0 , \sigma^{2}_{t}) \] \[ \mu_t = 0.001373 + -0.155505 \sigma^{2}_{t} \] \[ \sigma^{2}_{t} = 0.000018 + -0.062292 a^{2}_{t-1} + 0 \cdot \sigma^{2}_{t-1} + -0.062292 I(a^{2}_{t-1} < 0) a^{2}_{t-1} \]
報酬率是否存在風險溢酬(Risk Premium)與波動不對稱性(Asymmetric volatility effect)?
加入該模型之參數估計當作起始值,使rolling更快收斂。
ARMA_p = 0
ARMA_q = 0
GARCH_p = 0
GARCH_q = 1
# 取得係數作為起使值
coef_list = as.list(coef(modelfit2_2))
spec3 = ugarchspec(mean.model = list(armaOrder = c(ARMA_p,ARMA_q),
include.mean = TRUE,
archm = TRUE,
archpow = 1,
external.regressors = NULL,
archex = FALSE),
variance.model = list(model = "gjrGARCH",
garchOrder = c(GARCH_q, GARCH_p),
external.regressors = NULL,
variance.targeting = TRUE),
distribution.model = "norm",
start.pars = coef_list)
# rolling1 = ugarchroll(spec = spec3 ,
# data = Asset$Asset.LogReturn, # 放入全部的資料
# n.ahead = 1,
# forecast.length = OutSample,
# refit.every = 1,
# refit.window = "moving",
# solver = "hybrid",
# solver.control = solver_control,
# # fit.control = list(scale = 1),
# calculate.VaR = FALSE,
# parallel = TRUE,
# parallel.control = list(pkg = c("snowfall"), cores = 4),
# keep.coef = TRUE)
#
# rolling1 = resume(rolling1 , solver="gosolnp")
# save(rolling1 , file = 'C:/Users/amyhs/Desktop/碩士課程/時間序列/期中/rolling1.RData')
load('C:/Users/amyhs/Desktop/碩士課程/時間序列/期中/rolling1.RData')
convergence(rolling1)
## [1] 0
結果為0表示收斂。
coef(rolling1)[[1]]
## $index
## [1] "2017-12-29"
##
## $coef
## Estimate Std. Error t value Pr(>|t|)
## mu 1.366629e-03 0.001721836 0.7937045 0.4273675
## archm -1.555664e-01 0.413306718 -0.3763946 0.7066236
## alpha1 3.134366e-02 0.066370283 0.4722544 0.6367452
## gamma1 -6.229181e-02 0.069521147 -0.8960124 0.3702461
## omega 1.791456e-05 NA NA NA
由上表可知我們需要的estimated archm在[2,1]的位置,p-value of archm在[2,4]的位置;estimated gamma1在[4,1]的位置,p-value of gamma1在[4,4]的位置。
>> archm :表風險溢酬效應
>> gamma1:表波動不對稱效應
rolling1_time = as.Date(sapply(1:OutSample, function(x) coef(rolling1)[[x]]$index))
Asset$archm_estimate_rolling1 = NA
Asset$archm_pvalue_rolling1 = NA
Asset$archm_estimate_rolling1[WindowSize:(DataSize-1)] = sapply(1:OutSample, function(x) coef(rolling1)[[x]]$coef[2,1])
Asset$archm_pvalue_rolling1[WindowSize:(DataSize-1)] = sapply(1:OutSample, function(x) coef(rolling1)[[x]]$coef[2,4])
prem1 = cbind(Asset$archm_estimate_rolling1 , Asset$archm_pvalue_rolling1)
alpha = 0.05
prem1$h = alpha
prem1 = na.omit(prem1)
| 圖中可發現archm的p-value大多大於0.05,因此判斷參數估計不顯著,無風險溢酬的效果。 |
Asset$gamma_estimate_rolling1 = NA
Asset$gamma_pvalue_rolling1 = NA
Asset$gamma_estimate_rolling1[WindowSize:(DataSize-1)] = sapply(1:OutSample, function(x) coef(rolling1)[[x]]$coef[4,1])
Asset$gamma_pvalue_rolling1[WindowSize:(DataSize-1)] = sapply(1:OutSample, function(x) coef(rolling1)[[x]]$coef[4,4])
vol1 = cbind(Asset$gamma_estimate_rolling1 , Asset$gamma_pvalue_rolling1)
vol1$h = alpha
vol1 = na.omit(vol1)
| 圖中可看到僅有2018年2、3月gamma1的p-value < 0.05,有波動不對稱的現象;其餘時間gamma1的p-value > 0.05,便無明顯波動不對稱的現象。 |
S&P 500 指數日平均報酬率、S&P 500指數日報酬率平均真實波動 (realized volatility)、風險溢酬、波動不對稱性之關係
# 設定平均之天數
windowns = 20
# 計算平均真實波動
real_vol = function(sqare_r_t){
return(sqrt(sum(sqare_r_t)/windowns))
}
# S&P 500日平均報酬率
Asset$avg.return = rollapply(Asset$Asset.LogReturn , width = windowns , FUN = 'mean')
# 日報酬率平均真實波動 (realized volatility)
Asset$sqare_r = Asset$Asset.LogReturn^2
Asset$real.vol = rollapply(Asset$sqare_r , width = windowns , FUN = real_vol)
將畫圖資料合併
plot_data = cbind(Asset$archm_estimate_rolling1 , Asset$gamma_estimate_rolling1 , Asset$avg.return , Asset$real.vol)
plot_data = na.omit(plot_data)
嘗試用迴圈找出所有組合並畫圖,但跑不出來
# # 嘗試用迴圈找出所有組合,但做不出來
# set = t(combn(names(plot_data), 2))
#
# for (i in 1:6){
# data1 = plot_ly(data = as.data.frame(plot_data) , x = index(plot_data)) %>%
#
# add_lines(y = plot_data$set[,1][i] , type = "scatter" , mode = "lines" ,
# line = list(color = '#D79E4B') ,
# name = set[,1][i]) %>%
#
# layout(title = sprintf('%s v.s. %s' , set[,1][i] , set[,2][i]),
# xaxis = list(
# rangeselector = list(
# buttons = list(
# list(
# count = 3,
# label = "3m",
# step = "month",
# stepmode = "backward"),
# list(
# count = 6,
# label = "6m",
# step = "month",
# stepmode = "backward"),
#
# list(step = "all"))),
#
# rangeslider = list(type = "date")),
#
# yaxis = list(side = 'left' ,
# title = set[,1][i]),
# legend = list(x = 0., y = 0.55 , orientation = 'h'))
#
# # plot p-value
# data2 = plot_ly(data = as.data.frame(plot_data) , x = index(plot_data)) %>%
#
# add_lines(y = plot_data$set[,2][i] , type = "scatter" , mode = "lines" ,
# line = list(color = 'cornflowerblue') ,
# name = '') %>%
#
# layout(yaxis = list(side = 'left' ,
# title = set[,2][i]))
#
# fig_data = subplot(data1 , data2 , nrows = 2 , shareX = TRUE , margin = 0.05)
# offline(fig_data)
# break
# }
改成兩兩作圖
散佈圖
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
相關係數
cor(plot_data$avg.return , plot_data$real.vol)
## real.vol
## avg.return -0.739848
| 由散佈圖和相關係數(-0.74)來看,可以確定兩者為中度負相關。 |
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
相關係數
cor(plot_data$archm_estimate_rolling1 , plot_data$real.vol)
## real.vol
## archm_estimate_rolling1 -0.3158987
| 由散佈圖和相關係數(-0.32)來看,可以確定兩者為低度負相關。 |
。
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
相關係數
cor(plot_data$avg.return , plot_data$gamma_estimate_rolling1)
## gamma_estimate_rolling1
## avg.return -0.1488691
| 由散佈圖和相關係數(-0.15)來看,可以確定兩者為低度負相關。 |
。
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
相關係數
cor(plot_data$real.vol , plot_data$gamma_estimate_rolling1)
## gamma_estimate_rolling1
## real.vol -0.09533797
| 由散佈圖和相關係數(-0.10)來看,可以確定兩者為低度負相關。 |
。
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
相關係數
cor(plot_data$avg.return , plot_data$archm_estimate_rolling1)
## archm_estimate_rolling1
## avg.return 0.07314865
| 由散佈圖和相關係數(0.07)來看,可以確定兩者幾乎無相關。 |
。
加入星期一效應為解釋變數
由於解釋變數可以加在mean equation,也可以加在variable equation,因此以下假設三個模型做擬合,再找出最式模型。
>> 加入該模型之參數估計當作起始值,使rolling更快收斂
# 取得係數作為起使值
coef_list = as.list(coef(modelfit2_2))
spec4_1 = ugarchspec(mean.model = list(armaOrder = c(ARMA_p,ARMA_q),
include.mean = TRUE,
archm = TRUE,
archpow = 1,
external.regressors = Regressor_mean,
archex = FALSE),
variance.model=list(model = "gjrGARCH",
garchOrder = c(GARCH_q, GARCH_p),
external.regressors = NULL,
variance.targeting = TRUE),
distribution.model = "norm",
start.pars = coef_list)
modelfit4_1 = ugarchfit(spec = spec4_1,
data = Return_2017,
solver = "hybrid",
solver.control = solver_control
)
##
## Iter: 1 fn: -1019.5340 Pars: 0.0012924 -0.1627963 0.0005839 0.0313437 -0.0622919
## Iter: 2 fn: -1019.5340 Pars: 0.0012924 -0.1627963 0.0005839 0.0313437 -0.0622919
## solnp--> Completed in 2 iterations
modelfit4_1
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : gjrGARCH(1,0)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001292 0.001064 1.21494 0.22439
## archm -0.162796 0.243983 -0.66724 0.50462
## mxreg1 0.000584 0.000681 0.85709 0.39140
## alpha1 0.031344 0.088053 0.35596 0.72187
## gamma1 -0.062292 0.096120 -0.64806 0.51695
## omega 0.000018 NA NA NA
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001292 0.001521 0.84978 0.39545
## archm -0.162796 0.379274 -0.42923 0.66775
## mxreg1 0.000584 0.000734 0.79577 0.42617
## alpha1 0.031344 0.068203 0.45956 0.64583
## gamma1 -0.062292 0.069247 -0.89956 0.36835
## omega 0.000018 NA NA NA
##
## LogLikelihood : 1019.534
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.0839
## Bayes -8.0137
## Shibata -8.0847
## Hannan-Quinn -8.0557
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 4.637 0.03129
## Lag[2*(p+q)+(p+q)-1][2] 4.992 0.04087
## Lag[4*(p+q)+(p+q)-1][5] 5.614 0.11073
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.03736 0.8467
## Lag[2*(p+q)+(p+q)-1][2] 0.03952 0.9646
## Lag[4*(p+q)+(p+q)-1][5] 0.84340 0.8940
## d.o.f=1
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[2] 0.004253 0.500 2.000 0.9480
## ARCH Lag[4] 0.125808 1.397 1.611 0.9781
## ARCH Lag[6] 1.879858 2.222 1.500 0.7072
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 0.7503
## Individual Statistics:
## mu 0.01628
## archm 0.04579
## mxreg1 0.20027
## alpha1 0.34619
## gamma1 0.22136
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.9690 0.3335
## Negative Sign Bias 0.1987 0.8426
## Positive Sign Bias 0.2457 0.8061
## Joint Effect 2.2309 0.5259
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 44.22 0.0008806
## 2 30 60.35 0.0005579
## 3 40 59.92 0.0172242
## 4 50 79.08 0.0041612
##
##
## Elapsed time : 0.08563089
modelfit4_1_std_residual = modelfit4_1@fit$residuals/modelfit4_1@fit$sigma
par(mfcol = c(2,2) , mai=c(0.75,0.75,0.7,0.7))
acf(modelfit4_1_std_residual, main = sprintf('%s - Standardized Residual' , ticker))
pacf(modelfit4_1_std_residual , main = '')
acf((modelfit4_1_std_residual)^2. , main = sprintf('%s - Square of Standardized Residual' , ticker))
pacf((modelfit4_1_std_residual)^2. , main = '')
Weighted.Box.test(modelfit4_1_std_residual , lag = 10, type="Ljung-Box")
Weighted.Box.test(modelfit4_1_std_residual , lag = 20, type="Ljung-Box")
Weighted.Box.test(modelfit4_1_std_residual^2 , lag = 10, type="Ljung-Box")
Weighted.Box.test(modelfit4_1_std_residual^2 , lag = 20, type="Ljung-Box")
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_1_std_residual
## Weighted X-squared on Residuals for fitted ARMA process = 7.3121,
## Shape = 3.9286, Scale = 1.4000, p-value = 0.2244
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_1_std_residual
## Weighted X-squared on Residuals for fitted ARMA process = 12.173,
## Shape = 7.6829, Scale = 1.3667, p-value = 0.2948
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_1_std_residual^2
## Weighted X-squared on Residuals for fitted ARMA process = 3.0738,
## Shape = 3.9286, Scale = 1.4000, p-value = 0.8093
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_1_std_residual^2
## Weighted X-squared on Residuals for fitted ARMA process = 6.359,
## Shape = 7.6829, Scale = 1.3667, p-value = 0.8765
spec4_2 = ugarchspec(mean.model = list(armaOrder = c(ARMA_p,ARMA_q),
include.mean = TRUE,
archm = TRUE,
archpow = 1,
external.regressors = NULL,
archex = FALSE),
variance.model=list(model = "gjrGARCH",
garchOrder = c(GARCH_q, GARCH_p),
external.regressors = Regressor_var,
variance.targeting = TRUE),
distribution.model = "norm",
start.pars = coef_list)
modelfit4_2 = ugarchfit(spec = spec4_2,
data = Return_2017,
solver = "hybrid",
solver.control = solver_control
)
##
## Iter: 1 fn: -1019.1398 Pars: 0.00136295 -0.15562500 0.03134343 -0.06238857 0.00000001
## solnp--> Completed in 1 iterations
modelfit4_2
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : gjrGARCH(1,0)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001363 0.000761 1.791741 0.073174
## archm -0.155625 0.187709 -0.829078 0.407061
## alpha1 0.031343 0.036820 0.851253 0.394629
## gamma1 -0.062389 0.033335 -1.871558 0.061268
## vxreg1 0.000000 0.000000 0.047419 0.962179
## omega 0.000018 NA NA NA
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001363 0.001190 1.145780 0.251886
## archm -0.155625 0.297494 -0.523120 0.600891
## alpha1 0.031343 0.062740 0.499575 0.617375
## gamma1 -0.062389 0.034492 -1.808797 0.070483
## vxreg1 0.000000 0.000000 0.049586 0.960453
## omega 0.000018 NA NA NA
##
## LogLikelihood : 1019.14
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.0808
## Bayes -8.0106
## Shibata -8.0816
## Hannan-Quinn -8.0525
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 4.382 0.03631
## Lag[2*(p+q)+(p+q)-1][2] 4.842 0.04482
## Lag[4*(p+q)+(p+q)-1][5] 5.527 0.11596
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.02972 0.8631
## Lag[2*(p+q)+(p+q)-1][2] 0.05588 0.9518
## Lag[4*(p+q)+(p+q)-1][5] 0.96229 0.8682
## d.o.f=1
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[2] 0.05148 0.500 2.000 0.8205
## ARCH Lag[4] 0.23462 1.397 1.611 0.9497
## ARCH Lag[6] 2.07831 2.222 1.500 0.6627
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 0.9402
## Individual Statistics:
## mu 0.01518
## archm 0.04198
## alpha1 0.34038
## gamma1 0.21674
## vxreg1 0.52109
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.6758 0.4998
## Negative Sign Bias 0.0670 0.9466
## Positive Sign Bias 0.4028 0.6875
## Joint Effect 1.7206 0.6324
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 59.36 4.886e-06
## 2 30 73.50 9.913e-06
## 3 40 78.40 1.861e-04
## 4 50 93.82 1.223e-04
##
##
## Elapsed time : 0.02462983
modelfit4_2_std_residual = modelfit4_2@fit$residuals/modelfit4_2@fit$sigma
par(mfcol = c(2,2) , mai=c(0.75,0.75,0.7,0.7))
acf(modelfit4_2_std_residual, main = sprintf('%s - Standardized Residual' , ticker))
pacf(modelfit4_2_std_residual , main = '')
acf((modelfit4_2_std_residual)^2. , main = sprintf('%s - Square of Standardized Residual' , ticker))
pacf((modelfit4_2_std_residual)^2. , main = '')
Weighted.Box.test(modelfit4_2_std_residual , lag = 10, type="Ljung-Box")
Weighted.Box.test(modelfit4_2_std_residual , lag = 20, type="Ljung-Box")
Weighted.Box.test(modelfit4_2_std_residual^2 , lag = 10, type="Ljung-Box")
Weighted.Box.test(modelfit4_2_std_residual^2 , lag = 20, type="Ljung-Box")
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_2_std_residual
## Weighted X-squared on Residuals for fitted ARMA process = 7.1608,
## Shape = 3.9286, Scale = 1.4000, p-value = 0.2381
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_2_std_residual
## Weighted X-squared on Residuals for fitted ARMA process = 11.972,
## Shape = 7.6829, Scale = 1.3667, p-value = 0.3117
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_2_std_residual^2
## Weighted X-squared on Residuals for fitted ARMA process = 3.292,
## Shape = 3.9286, Scale = 1.4000, p-value = 0.7769
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_2_std_residual^2
## Weighted X-squared on Residuals for fitted ARMA process = 6.6588,
## Shape = 7.6829, Scale = 1.3667, p-value = 0.8528
spec4_3 = ugarchspec(mean.model = list(armaOrder = c(ARMA_p,ARMA_q),
include.mean = TRUE,
archm = TRUE,
archpow = 1,
external.regressors = Regressor_mean,
archex = FALSE),
variance.model=list(model = "gjrGARCH",
garchOrder = c(GARCH_q, GARCH_p),
external.regressors = Regressor_var,
variance.targeting = TRUE),
distribution.model = "norm",
start.pars = coef_list)
modelfit4_3 = ugarchfit(spec = spec4_3,
data = Return_2017,
solver = "hybrid",
solver.control = solver_control
)
##
## Iter: 1 fn: -1019.5304 Pars: 0.00130191 -0.16192470 0.00060575 0.03134342 -0.06240600 0.00000001
## Iter: 2 fn: -1019.5371 Pars: 0.00129147 -0.16253869 0.00055478 0.03134150 -0.06321167 0.00000001
## solnp--> Completed in 2 iterations
modelfit4_3
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : gjrGARCH(1,0)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001291 0.000749 1.725177 0.084496
## archm -0.162539 0.182684 -0.889724 0.373614
## mxreg1 0.000555 0.000682 0.814025 0.415630
## alpha1 0.031342 0.035971 0.871302 0.383589
## gamma1 -0.063212 0.031658 -1.996709 0.045857
## vxreg1 0.000000 0.000000 0.052891 0.957819
## omega 0.000018 NA NA NA
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.001291 0.001063 1.214491 0.224560
## archm -0.162539 0.282457 -0.575446 0.564989
## mxreg1 0.000555 0.000738 0.751823 0.452158
## alpha1 0.031342 0.060056 0.521875 0.601758
## gamma1 -0.063212 0.036135 -1.749315 0.080237
## vxreg1 0.000000 0.000000 0.058339 0.953479
## omega 0.000018 NA NA NA
##
## LogLikelihood : 1019.537
##
## Information Criteria
## ------------------------------------
##
## Akaike -8.0760
## Bayes -7.9917
## Shibata -8.0771
## Hannan-Quinn -8.0421
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 4.638 0.03127
## Lag[2*(p+q)+(p+q)-1][2] 4.998 0.04071
## Lag[4*(p+q)+(p+q)-1][5] 5.620 0.11036
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.02979 0.8630
## Lag[2*(p+q)+(p+q)-1][2] 0.03241 0.9703
## Lag[4*(p+q)+(p+q)-1][5] 0.83864 0.8950
## d.o.f=1
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[2] 0.005164 0.500 2.000 0.9427
## ARCH Lag[4] 0.129152 1.397 1.611 0.9773
## ARCH Lag[6] 1.882594 2.222 1.500 0.7066
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.3929
## Individual Statistics:
## mu 0.01565
## archm 0.04309
## mxreg1 0.19854
## alpha1 0.29907
## gamma1 0.17723
## vxreg1 0.68141
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.49 1.68 2.12
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.9616 0.3372
## Negative Sign Bias 0.1786 0.8584
## Positive Sign Bias 0.2486 0.8039
## Joint Effect 2.2375 0.5246
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 44.06 0.0009266
## 2 30 62.03 0.0003432
## 3 40 60.87 0.0140467
## 4 50 83.46 0.0015555
##
##
## Elapsed time : 0.06516385
modelfit4_3_std_residual = modelfit4_3@fit$residuals/modelfit4_3@fit$sigma
par(mfcol = c(2,2) , mai=c(0.75,0.75,0.7,0.7))
acf(modelfit4_3_std_residual , main = sprintf('%s - Standardized Residual' , ticker))
pacf(modelfit4_3_std_residual , main = '')
acf((modelfit4_3_std_residual)^2. , main = sprintf('%s - Square of Standardized Residual' , ticker))
pacf((modelfit4_3_std_residual)^2. , main = '')
Weighted.Box.test(modelfit4_3_std_residual , lag = 10, type="Ljung-Box")
Weighted.Box.test(modelfit4_3_std_residual , lag = 20, type="Ljung-Box")
Weighted.Box.test(modelfit4_3_std_residual^2 , lag = 10, type="Ljung-Box")
Weighted.Box.test(modelfit4_3_std_residual^2 , lag = 20, type="Ljung-Box")
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_3_std_residual
## Weighted X-squared on Residuals for fitted ARMA process = 7.3129,
## Shape = 3.9286, Scale = 1.4000, p-value = 0.2243
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_3_std_residual
## Weighted X-squared on Residuals for fitted ARMA process = 12.174,
## Shape = 7.6829, Scale = 1.3667, p-value = 0.2948
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_3_std_residual^2
## Weighted X-squared on Residuals for fitted ARMA process = 3.0708,
## Shape = 3.9286, Scale = 1.4000, p-value = 0.8098
##
##
## Weighted Ljung-Box test (Gamma Approximation)
##
## data: modelfit4_3_std_residual^2
## Weighted X-squared on Residuals for fitted ARMA process = 6.3642,
## Shape = 7.6829, Scale = 1.3667, p-value = 0.8762
| 模 型 | ACF/PACF | Weighted Ljuang-box Test | BIC | 是否選用 |
|---|---|---|---|---|
| 模型一 | 無自我相關 | 通過 | -8.0137 | O |
| 模型二 | 無自我相關 | 通過 | -8.0106 | X |
| 模型三 | 無自我相關 | 通過 | -7.9917 | X |
| 三個模型之殘差皆可視為white noise,但最後選擇BIC最小之模型為最適模型。 |
Rolling Analysis
# rolling2 = ugarchroll(spec = spec4_1 ,
# data = Asset$Asset.LogReturn, # 放入全部的資料
# n.ahead = 1,
# forecast.length = OutSample,
# refit.every = 1,
# refit.window = "moving",
# solver = "hybrid",
# solver.control = solver_control,
# fit.control = list(scale = 1),
# calculate.VaR = FALSE,
# parallel = TRUE,
# parallel.control = list(pkg = c("snowfall"), cores = 4),
# keep.coef = TRUE)
#
# rolling2 = resume(rolling2 , solver="gosolnp")
# save(rolling2 , file = 'C:/Users/amyhs/Desktop/碩士課程/時間序列/期中/rolling2.RData')
load('C:/Users/amyhs/Desktop/碩士課程/時間序列/期中/rolling2.RData')
convergence(rolling2)
## [1] 0
結果為0表示收斂。
coef(rolling2)[[1]]
## $index
## [1] "2017-12-29"
##
## $coef
## Estimate Std. Error t value Pr(>|t|)
## mu 1.292397e-03 0.0015208635 0.8497785 0.3954482
## archm -1.627963e-01 0.3792742448 -0.4292312 0.6677550
## mxreg1 5.838820e-04 0.0007337308 0.7957715 0.4261649
## alpha1 3.134368e-02 0.0682032689 0.4595627 0.6458301
## gamma1 -6.229185e-02 0.0692469409 -0.8995611 0.3683539
## omega 1.790323e-05 NA NA NA
由上表可知我們需要的estimated archm在[2,1]的位置,p-value of archm在[2,4]的位置;estimated gamma1在[5,1]的位置,p-value of gamma1在[5,4]的位置。
rolling2_time = as.Date(sapply(1:OutSample, function(x) coef(rolling2)[[x]]$index))
Asset$archm_estimate_rolling2 = NA
Asset$archm_pvalue_rolling2 = NA
Asset$archm_estimate_rolling2[WindowSize:(DataSize-1)] = sapply(1:OutSample, function(x) coef(rolling2)[[x]]$coef[2,1])
Asset$archm_pvalue_rolling2[WindowSize:(DataSize-1)] = sapply(1:OutSample, function(x) coef(rolling2)[[x]]$coef[2,4])
prem2 = cbind(Asset$archm_estimate_rolling2 , Asset$archm_pvalue_rolling2)
alpha = 0.05
prem2$h = alpha
prem2 = na.omit(prem2)
# plot estimated
p5 = plot_ly(data = as.data.frame(prem2) , x = index(prem2)) %>%
add_lines(y = ~archm_estimate_rolling2 , type = "scatter" , mode = "lines" ,
line = list(color = '#D79E4B') ,
name = 'Estimated archm (rolling2)') %>%
layout(title = sprintf('%s ARMA-GJR-Garch:archm (with regressor)' , ticker),
xaxis = list(
rangeselector = list(
buttons = list(
list(
count = 3,
label = "3m",
step = "month",
stepmode = "backward"),
list(
count = 6,
label = "6m",
step = "month",
stepmode = "backward"),
list(step = "all"))),
rangeslider = list(type = "date")),
yaxis = list(side = 'left' ,
title = 'Estimated archm'),
legend = list(x = 0., y = 0.55 , orientation = 'h'))
# plot p-value
p6 = plot_ly(data = as.data.frame(prem2) , x = index(prem2)) %>%
add_lines(y = ~archm_pvalue_rolling2 , type = "scatter" , mode = "lines" ,
line = list(color = 'cornflowerblue') ,
name = 'p-value of archm') %>%
add_lines(y = ~h , type = "scatter" , mode = "lines" ,
line = list(color = 'red') , name = sprintf("alpha(%.2f)" , alpha))
P3 = subplot(p5 , p6 , nrows = 2 , shareX = TRUE , titleY = TRUE , margin = 0.08)
offline(P3)
# htmlwidgets::saveWidget( P3 , 'C:/Users/amyhs/Desktop/碩士課程/時間序列/模擬/Figure/Risk Premium (rolling2).html')
| 圖中可發現archm的p-value大多大於0.05,因此判斷參數估計不顯著,無風險溢酬的效果。 |
Asset$gamma_estimate_rolling2 = NA
Asset$gamma_pvalue_rolling2 = NA
Asset$gamma_estimate_rolling2[WindowSize:(DataSize-1)] = sapply(1:OutSample, function(x) coef(rolling2)[[x]]$coef[5,1])
Asset$gamma_pvalue_rolling2[WindowSize:(DataSize-1)] = sapply(1:OutSample, function(x) coef(rolling2)[[x]]$coef[5,4])
vol2 = cbind(Asset$gamma_estimate_rolling2 , Asset$gamma_pvalue_rolling2)
vol2$h = alpha
vol2 = na.omit(vol2)
# plot estimated
p7 = plot_ly(data = as.data.frame(vol2) , x = index(vol2)) %>%
add_lines(y = ~gamma_estimate_rolling2 , type = "scatter" , mode = "lines" ,
line = list(color = '#D79E4B') ,
name = 'Estimated Gamma') %>%
layout(title = sprintf('%s ARMA-GJR-Garch:Gamma (with regressor)' , ticker),
xaxis = list(
rangeselector = list(
buttons = list(
list(
count = 3,
label = "3m",
step = "month",
stepmode = "backward"),
list(
count = 6,
label = "6m",
step = "month",
stepmode = "backward"),
list(step = "all"))),
rangeslider = list(type = "date")),
yaxis = list(side = 'left' ,
title = 'Estimated Gamma'),
legend = list(x = 0., y = 0.55 , orientation = 'h'))
# plot p-value
p8 = plot_ly(data = as.data.frame(vol2) , x = index(vol2)) %>%
add_lines(y = ~gamma_pvalue_rolling2 , type = "scatter" , mode = "lines" ,
line = list(color = 'cornflowerblue') ,
name = 'p-value of Gamma') %>%
add_lines(y = ~h , type = "scatter" , mode = "lines" ,
line = list(color = 'red') , name = sprintf("alpha(%.2f)" , alpha))
P4 = subplot(p7 , p8 , nrows = 2 , shareX = TRUE , titleY = TRUE , margin = 0.08)
offline(P4)
# htmlwidgets::saveWidget( P4 , 'C:/Users/amyhs/Desktop/碩士課程/時間序列/模擬/Figure/Asymmetric Parameter (rolling2).html')
| 圖中可看到2018年7月到11月gamma的p-value幾乎大於0.05,無顯著波動不對稱的現象;但其餘時間gamma的p-value幾乎大於0.05,有顯著波動不對稱的現象。 |
dataset = na.omit(as.data.frame(Asset))
p9 = plot_ly(data = dataset , x = index(dataset)) %>%
add_lines(y = ~archm_estimate_rolling1 , type = "scatter" , mode = "lines" ,
line = list(color = '#D79E4B') ,
name = 'without regressor') %>%
add_lines(y = ~archm_estimate_rolling2 , type = "scatter" , mode = "lines" ,
line = list(color = 'cornflowerblue') ,
name = 'with regressor' , yaxis = "y2") %>%
layout(title = sprintf('%s Archm (without regressor vs with regressor)' , ticker),
xaxis = list(
rangeselector = list(
buttons = list(
list(
count = 3,
label = "3m",
step = "month",
stepmode = "backward"),
list(
count = 6,
label = "6m",
step = "month",
stepmode = "backward"),
list(step = "all"))),
rangeslider = list(type = "date")),
yaxis = list(side = 'left' ,
title = 'without regressor'),
yaxis2 = list(title = 'with regressor' ,
overlaying = "y", side = "right"),
legend = list(x = 0.55, y = 0.15 , orientation = 'h'))
offline(p9)
p10 = plot_ly(data = dataset , x = index(dataset)) %>%
add_lines(y = ~gamma_estimate_rolling1 , type = "scatter" , mode = "lines" ,
line = list(color = '#D79E4B') ,
name = 'without regressor') %>%
add_lines(y = ~gamma_estimate_rolling2 , type = "scatter" , mode = "lines" ,
line = list(color = 'cornflowerblue') ,
name = 'with regressor' , yaxis = "y2") %>%
layout(title = sprintf('%s Gamma (without regressor vs with regressor)' , ticker),
xaxis = list(
rangeselector = list(
buttons = list(
list(
count = 3,
label = "3m",
step = "month",
stepmode = "backward"),
list(
count = 6,
label = "6m",
step = "month",
stepmode = "backward"),
list(step = "all"))),
rangeslider = list(type = "date")),
yaxis = list(side = 'left' ,
title = 'without regressor'),
yaxis2 = list(title = 'with regressor' ,
overlaying = "y", side = "right"),
legend = list(x = 0.55, y = 0.15 , orientation = 'h'))
offline(p10)
| 參數 | 無解釋變數 | 有解釋變數 |
|---|---|---|
| archm | 不顯著 | 不顯著 |
| gamma | 不顯著 | 顯著 |
加入解釋變數後,gamma效應大部分時間都顯著,因此判斷星期一效應加入模型有助模型的解釋能力。